from pathlib import Path

import matplotlib.pyplot as plt
import numpy as np
import scipy.integrate as integrate
from docx.shared import Cm
from docxtpl import DocxTemplate, InlineImage
from scipy.misc import derivative

OUTPUT_FOLDER = Path("dz/heat_resistance")
DOCX_TEMPLATE_PATH = OUTPUT_FOLDER.joinpath("docx_template.docx")
REPORT_OUTPUT_PATH = OUTPUT_FOLDER.joinpath("report.docx")
TEMPORARY_IMAGE_PATH_1 = OUTPUT_FOLDER.joinpath("image1.png")
TEMPORARY_IMAGE_PATH_2 = OUTPUT_FOLDER.joinpath("image2.png")
TEMPORARY_IMAGE_PATH_3 = OUTPUT_FOLDER.joinpath("image3.png")
TEMPORARY_IMAGE_PATH_4 = OUTPUT_FOLDER.joinpath("image4.png")
TEMPORARY_IMAGE_PATH_5 = OUTPUT_FOLDER.joinpath("image5.png")
TEMPORARY_IMAGE_PATH_6 = OUTPUT_FOLDER.joinpath("image6.png")
# Исходные данные
E = 6  # Модуль Юнга, МПа
a = 0.3  # Внутренний радиус, м
b = 1.5  # Внешний радиус, м
v = 0.4  # Коэффициент Пуассона

Ta = 300  # Температура внутренней стенки, °С
Tb = 100  # Температура внешней стенки, °С
al_t = 10 ** (-5)  # Коэффициент теплопередачи
Bia = 1 / a  # Число Био для внутренней стенки, 1/м
Bib = 1 / b  #  Число Био для внешней стенки, 1/м
m = E / (2 * (1 + v))  # Параметр Ламе
la = v * E / ((1 + v) * (1 - 2 * v))  # Параметр Ламе
D2 = 1  # Тензор скорости деформаций


def f2(dmnless_time: float) -> float:
    """Часть функции напряжений

    Параметры
    ----------
    dmnless_time : float
        Время, сек

    Возвращает
    -------
    float
        Результат
    """
    result = (
        2
        - (1 + v) * np.exp(-(1 - v) * dmnless_time / (1 + v))
        - (1 - v) * np.exp(-dmnless_time)
    )
    return result


def g_mp_r(r: float) -> float:
    """Часть функции температурных деформаций, помноженная на r

    Параметры
    ----------
    r : float
        Радиус, м

    Возвращает
    -------
    float
        Реузльтат
    """
    result = al_t * (c1 * np.log10(r) + c2)
    return result * r


def Q(r: float) -> float:
    """Часть функции напряжения

    Параметры
    ----------
    r : float
        Радиус, м

    Возвращает
    -------
    float
        Результат
    """
    result = integrate.quad(g_mp_r, a, r)[0]
    return result


def sig_rr(r: float, dmnless_time: float) -> float:
    """Функция радиальных напряжений

    Параметры
    ----------
    r : float
        Радиус, м
    dmnless_time : float
        Время, с

    Возвращает
    -------
    float
        Напряжение, МПа
    """
    result = (
        f2(dmnless_time)
        * (m * (3 * la + 2 * m) / (la + 2 * m))
        * ((1 - a**2 / r**2) / (b**2 - a**2) * Q(b) - Q(r) / r**2)
    )
    return result


def sig_phi_phi(r: float, dmnless_time: float) -> float:
    """Функция тангенциальных напряжений

    Параметры
    ----------
    r : float
        Радиус, м
    dmnless_time : float
        Время, сек

    Возвращает
    -------
    float
        Напряжение, МПа
    """
    result = (
        f2(dmnless_time)
        * (m * (3 * la + 2 * m) / (la + 2 * m))
        * (
            (1 - a**2 / r**2) / (b**2 - a**2) * Q(b)
            + Q(r) / r**2
            - g_mp_r(r) / r
        )
    )
    return result


def save_plot(
    x_list: list,
    y_list: list,
    title: str,
    x_label: str,
    y_label: str,
    path: any,
    isOneScaled: bool = False,
) -> None:
    """Функция сохранения графика

    Параметры
    ----------
    x_list : list
        Список x составляющей точек
    y_list : list
        Список y составляющей точек
    title : str
        Заголовок графика
    x_label : str
        Заголовок оси X
    y_label : str
        Заголовок оси Y
    path : any
        Путь сохранения графика
    isOneScaled : bool, optional
        Нужно ли масштабировать график по оси Y от -1 до 1(для малых величин)
    """
    fig = plt.figure()
    if isOneScaled:
        plt.ylim([-1, 1])
    plt.grid(True)
    plt.title(title)
    plt.xlabel(x_label)
    plt.ylabel(y_label)
    plt.plot(x_list, y_list, color="red")
    fig.savefig(path, dpi=200)


def deriv_sig_zz(r: float, dmnless_time: float) -> float:
    """Функция производной осевых напряжений

    Параметры
    ----------
    r : float
        Радиус, м
    dmnless_time : float
        Время, сек

    Возвращает
    -------
    float
        Результат
    """
    result = -derivative(sig_rr, 1, dx=1e-6, args=(dmnless_time,)) - derivative(
        sig_phi_phi, 1, dx=1e-6, args=(dmnless_time,)
    )
    return result


def sig_phi_phi_2(dmnless_time: float, r: float) -> float:
    """Функция тангенциальных напряжений

    Параметры
    ----------
    dmnless_time : float
        Время, сек
    r : float
        Радиус, м

    Возвращает
    -------
    float
        Напряжение, МПа
    """
    result = (
        f2(dmnless_time)
        * (m * (3 * la + 2 * m) / (la + 2 * m))
        * (
            (1 - a**2 / r**2) / (b**2 - a**2) * Q(b)
            + Q(r) / r**2
            - g_mp_r(r) / r
        )
    )
    return result


def sig_zz(r: float, dmnless_time: float) -> float:
    """Функция осевых напряжений

    Параметры
    ----------
    r : float
        Радиус, м
    dmnless_time : float
        Время, сек

    Возвращает
    -------
    float
        Напряжение, МПа
    """
    result = (
        (
            la
            * (sig_rr(r, dmnless_time) + sig_phi_phi(r, dmnless_time))
            / (2 * m * (3 * la + 2 * m))
            - g_mp_r(r) / r
            - 2 * D2 * deriv_sig_zz(r, dmnless_time)
        )
        * m
        * (3 * la + 2 * m)
        / (la + m)
    )
    return result


def e_phi_phi(r: float, dmnless_time: float) -> float:
    """Функция тангенциальных деформаций

    Parameters
    ----------
    r : float
        Радиус, м
    dmnless_time : float
        Время, сек

    Returns
    -------
    float
        Деформация
    """
    result = (
        (
            sig_phi_phi(r, dmnless_time)
            - la
            / (3 * la + 2 * m)
            * (
                sig_phi_phi(r, dmnless_time)
                + sig_rr(r, dmnless_time)
                + sig_zz(r, dmnless_time)
            )
        )
        / (2 * m)
        + 2 * D2 * derivative(sig_phi_phi_2, 1, dx=1e-6, args=(r,))
        + g_mp_r(r) / r
    )
    return result * r


# Первая часть. Построение графика температурного поля
# -C1 == Bia * (Ta - C1*log(a)-C2)
# -C1 == Bib * (C1 * Log(b) + C2 - Tb)

A = np.array([[Bia * np.log10(a) - 1, Bia], [1 + Bib * np.log10(b), Bib]])

B = np.array([Bia * Ta, Bib * Tb])

c1, c2 = np.linalg.solve(A, B)
x_list = []
y_list = []
step = 0.01
r = a
end_r = b
while r <= end_r:
    x_list.append(r)
    y_list.append(c1 * np.log10(r) + c2)
    r += step

save_plot(
    x_list,
    y_list,
    "Температурное поле",
    "Радиус, м",
    "Температура, °С",
    TEMPORARY_IMAGE_PATH_1,
)

# Вторая часть. Распределение напряжений.
x_list = []
y_list = []
step = 0.01
r = a
end_r = b
while r <= end_r:
    x_list.append(r)
    y_list.append(sig_rr(r, 0))
    r += step

save_plot(
    x_list,
    y_list,
    "σrr (t->0)",
    "Радиус, м",
    "Напряжение, МПа",
    TEMPORARY_IMAGE_PATH_2,
    True,
)

x_list = []
y_list = []
step = 0.01
r = a
end_r = b
while r <= end_r:
    x_list.append(r)
    y_list.append(sig_rr(r, np.inf))
    r += step

save_plot(
    x_list, y_list, "σrr (t->∞)", "Радиус, м", "Напряжение, МПа", TEMPORARY_IMAGE_PATH_3
)

x_list = []
y_list = []
step = 0.01
r = a
end_r = b
while r <= end_r:
    x_list.append(r)
    y_list.append(sig_phi_phi(r, 0))
    r += step

save_plot(
    x_list,
    y_list,
    "σφφ (t->0)",
    "Радиус, м",
    "Напряжение, МПа",
    TEMPORARY_IMAGE_PATH_4,
    True,
)

x_list = []
y_list = []
step = 0.01
r = a
end_r = b
while r <= end_r:
    x_list.append(r)
    y_list.append(sig_phi_phi(r, np.inf))
    r += step

save_plot(
    x_list, y_list, "σφφ (t->∞)", "Радиус, м", "Напряжение, МПа", TEMPORARY_IMAGE_PATH_5
)

# Третья часть. Распределение радиальных перемещений
x_list = []
y_list = []
step = 0.01
r = a
end_r = b
while r <= end_r:
    x_list.append(r)
    y_list.append(e_phi_phi(r, 0))
    r += step

save_plot(
    x_list, y_list, "urr (t->0)", "Радиус, м", "Перемещение, м", TEMPORARY_IMAGE_PATH_6
)


docx_template = DocxTemplate(DOCX_TEMPLATE_PATH)
image_1 = InlineImage(docx_template, str(TEMPORARY_IMAGE_PATH_1), Cm(12))
image_2 = InlineImage(docx_template, str(TEMPORARY_IMAGE_PATH_2), Cm(12))
image_3 = InlineImage(docx_template, str(TEMPORARY_IMAGE_PATH_3), Cm(12))
image_4 = InlineImage(docx_template, str(TEMPORARY_IMAGE_PATH_4), Cm(12))
image_5 = InlineImage(docx_template, str(TEMPORARY_IMAGE_PATH_5), Cm(12))
image_6 = InlineImage(docx_template, str(TEMPORARY_IMAGE_PATH_6), Cm(12))
context = {
    "E": str(E),
    "a": str(a),
    "b": str(b),
    "v": str(v),
    "Ta": str(Ta),
    "Tb": str(Tb),
    "al_t": str(al_t),
    "image1": image_1,
    "image2": image_2,
    "image3": image_3,
    "image4": image_4,
    "image5": image_5,
    "image6": image_6,
}
docx_template.render(context)
docx_template.save(REPORT_OUTPUT_PATH)
